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Abstract. In this paper, we study a numerical approximation for a class of 
stationary states for reaction-diffusion system with m densities having disjoint 
support, which are governed by a minimization problem. We use quantitative 
properties of both solutions and free boundaries to derive our scheme. Fur- 
thermore, the proof of convergence of the numerical method is given in some 
particular cases. We also apply our numerical simulations for the spatial segre- 
gation limit of diffusive Lotka-Volterra models in presence of high competition 
and inhomogcncous Dirichlet boundary conditions. We discuss numerical im- 
plementations of the resulting approach and present computational tests. 



1. Introduction 

In recent years there have been intense studies of spatial segregation for 
reaction-diffusion systems . The existence of spatiaUy inhomogeneous solu- 
tions for competition models of Lotka-Volterra type in the case of two and 
more competing densities have been considered |CTV061 ICTV061 ICDHOTa] 
IDan951 IDZ02| Squ08| . The objective in this paper is to study numerical 



solutions for two classes of possible segregation states. The first class is re- 
lated with an arbitrary number of competing densities, which are governed 
by a minimization problem. 

Let ri C M", (n > 2) be a connected and bounded domain with smooth 
boundary and m be a fixed integer. We consider the steady-states of m 
competing species coexisting in the same area Let Ui{x) denotes the pop- 
ulation density of the i^^ component with the internal dynamic prescribed 
by fi{x). Here we assume that fi is uniformly continuous and fi{x) > 0. 

We call the m-tuple U = (ui, • • • , Um) G {H^{Q,))^, segregated states if 

Ui{x) • Uj{x) = 0, a.e. for i j, x £ Q. 
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Problem (A): Consider the following minimization problem 
(1.1) Minimize S(mi, ••• ,tim) = j ^ QlViijp + /jUj^ dx, 
over the set 

S = {{ui, . . . ,Um) € {H^i^))"" ■ Ui >0,Ui ■ Uj = 0,Ui = (pi on d^l}. 

Here (pi G H2 (c^O) with property (pi ■ (pj = 0, (/>j > on the boundary 90. 
Also we assume that fi is uniformly continuous and fi{x) > 0. 

Problem (B): Our second problem which appears in the study of pop- 
ulation ecology is the case when high competitive interactions between dif- 
ferent species occurs. As the rate of interaction of two different species 
goes to infinity, then competition-diffusion systems shows a limiting config- 
uration w ith segregated state. We refer the reader to |CTV05bl [CDHOTbt 
ICDH+n4[ IDD941 IDan95[ IDHMP991 IDZn2j and in particular to |DHMP99j 
for models involving Dirichlet boundary data. A complete analysis of the 
stationary case has been studied in |CTV05bJ . Also numerical simulation 
for the spatial segregation limit of two diffusive Lotka-Volterra models in 
presence of strong competition and inhomogeneous Dirichlet boundary con- 
ditions is provided in |SZ08] . In the work of [SZ08J the authors solve problem 
for small s and then let e — > 0, while in our work we use the qualitative 
properties of the limiting problem. Also unlike the results in |SZ08| . where 
the authors provided only simulations of their proposed algorithm, we give 
numerical consistent variational system with strong interaction, and provide 
disjointness condition of populations during iterations of the scheme. More- 
over by discussing those two problems we show that the proposed idea can 
be generalized for two or more species that competing each other. 

Let di, A be positive numbers. Consider the following system of m differ- 
ential equations 

\ Ui{x,y) = (pi{x,y) ondil, 

for i = 1, • • • , m, where again (pi G 2 (90) with property (pi ■ (pj = 0, (pi > 
on the boundary dO,. Our aim is to present numerical approximation for 
this system as e — )■ 0. This system can be viewed as a steady state of the 
following system in the case that boundary values are time independent. 



{^Ui - diAui = Xui{l - Ui) - luiY.jjkiu'] in X (0,oo), 
Ui{x, y, t) = (pi{x, y, t) on dQ x (0, 00), 

Ui{x,y,0) = Ui^o{x,y) in 0, 

for i = 1, • • • ,m. 

One of the interesting results which relates these two problems is given in 
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[CTVOSaj . Consider the following reaction-diffusion system of three com- 
peting species: 

(1.4) Aui = -Ui Uj, Ui > 0, in Ui = (pi, on dil. i = 1,2, 3, 

where we have the same assumptions on the boundary values In |CTV05a] 
it is shown the uniqueness of the limiting configuration as e — )• on a planar 
domain, with appropriate boundary conditions. Moreover, it was shown that 
the corresponding minimization problem admits a unique solution. Further- 
more, the limiting configuration minimizes the following energy 

^=l ^ 

over the set 5 = {ui E H^{Q) : Ui > 0,Ui ■ uj = 0, Uj = (pi ondCl, i = 1,2, 3}. 



To see the numerical approximations for the system of equations of (1.4) we 
refer to |Boz09j . 

2. Basic facts for Problem (A) 



In this section we will see that the solution for problem (1.1) satisfies a free 
boundary problem. In order to prove the existence of the minimizer we are 
going to apply the following classical theorem due to |Str90] . 

Theorem 2.1. Let V be a reflexive Banach space with norm \\ ■ \\, and 
M C V be a weakly closed subset of V. Suppose E : M ^ R is coercive on 
M with respect to V, that is 

i) E{u) — 7- oo as \\u\\ oo, u £ M and E is weakly lower semi- 
continuous on M with respect to V, that is 

ii) for any u G M, any sequence (um) in M such that Um u weakly 
in V there holds E{u) < liminf£'(um). 

Then E is bounded from below on M and attains its minimum in M. 

The main existence and uniqueness result for Problem (A) reads as fol- 
lows. 

Proposition 2.2. Under the assumptions in Problem (A), there exist a 



minimizer to (1.1), and it is unique. 



Proof. It is easy to see that the minimizer E{ui, ■ ■ ■ , Um) is coercive over 
the closed set S, and lower semi-continuous on S with respect to the space 



{H^(^l))"^. Thus the existence follows directly from theorem 2.1 stated above. 
For the proof of uniqueness we are using the same arguments as in [CTVOSbl 
Theorem 4.1]. Suppose that there exist two different minimizers U = 
(ui, • • • , Um) and V = {vi, ■ ■ ■ , Vm) of ( |1.1[ ) such that 

(2.1) E{ui, ■■■ , Um) = E{vi, ■■■ ,Vm) = C. 
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Define new functions Ui ^ Vi 9jS 

Ui{x) = Ui{x) -^Uk{x), 

Vi{x) = Vi{x) - "^Vkix), 

and let 

Wi{x) = -max {ui{x) +Vi{x),0) . 

Also 0.i denotes a positivity set of Wi, that is 

Qi = {x eU: Wi{x) > 0}. 

It is easy to see that Wi > 0, Wi ■ Wj = for i ^ j and Wi = (pi on and 
further we have 

E{wi,--- ,Wm)= / "^{-{Vwil^ + fiWi)dx. 

•^^ i=l 

Using that Ui ■ Uj = and Vi ■ Vj = for i ^ j, we can obtain the following 
estimate: 

/ ^2^^'^'^'^dx = Y^ / -\Vui + Vvifdx < / -{\Vuif + \Vvi\'^) dx 



< 



The potential part can also be estimated as follows 

/ ^fiWidx= j ^-fim&yi{ui{x)+Vi{x),Q)dx < j ^-fi{ui+Vi) dx, 

where in the last inequality wc have Tiscd the fact that fi is positive, {i = 
1, • • • ,m). Finally, by adding together two last inequalities we obtain 

E{wi,--- ,Wm) < ^[E{ui,--- ,Um) + E{vi, - ■ ■ ,Vm)] = C, 

which is a contradiction. □ 

In this part we state some results that will be used in the construction of 
our numerical scheme. These results show that the minimizer of the varia- 
tional problem satisfies a differential inequality, stated in the next lemma. 

Lemma 2.3. Let U = {ui,--- ,Um) be a minimizer of Problem (A), then 
the following holds in the sense of distributions. 

^Ui > fiix)X{ui>0}- 
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Proof. We should show that for each i = l,---m, and each test function 
(f) G C~(0) we have 

/ Vui ■ + fiX{ui>o}(l>dx < 0. 
Jn 

For < £ << 1, and fixed i we define new functions (vi, • • • Vm) as follows: 

Vi = Ui, V2 = U2,...,Vi = {Ui - ecf))'^, ...,Vm= Um- 

It is easy to see that, 

Vi ■ Vj = 0, whenever i ^ j, and Vi = (f)i on the boundary of il. 
Denote V = {vi, ■ ■ ■ Vm). We have E{U) < E{V), therefore 

< E{V) - E{U) = I \V{ui - scP)+\^ - \Vui\^ + fiiim - £0)+ - Ui)cl)dx 
Jn 

= / |V(?Xi -£^)+p - |VMipdx+ / fiX{ui>o}i{ui- e(py -Ui)4)dx 
Jn Jn 

X{ui>o}4') dx + o(£). 

Jn 



< 
Thus 



/ Vui ■ V4> + fiX{ui>o}<t> dx < 0. 
Jn 

□ 

For any minimizer U = (ui, ■ ■ ■ define the positive set of ui as 

Q.i = {xe^: Ui{x) > 0} . 

Definition 2.4. The multiplicity of a point a; G O is 

m{x) = card {i : measiyii fl B{x, r)) > for some r > 0} . 

The interface between two densities is defined as 

Fjj = d^i n d^j n{xe^: m{x) = 2}. 

Our numerical scheme is based on the following properties, which are 
straightforward to verify. 

Corollary 2.5. Assume that xq € ^ then the following holds: 

1) If m{xo) = 0, then there is r > such that for every i = 1, - ■ ■ m; Uj = 

on B{xo, r). 

2) If m{xQ) = 1, then there are i and r > such that in B{xo,r) 

Aui = fi, Uj = forj^i. 

3) If m{xo) = 2, then there are i,j and r > such that for every k and 
k ^ i,j we have Uk = and 

A{ui - Uj) = fixni - fjXnj in B{xo, r). 
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2.1. Special cases of Problem (A) 

We note that the one phase obstacle problem and the two-phase obstacle 
problem are special cases of Problem (A) for m = 1 and m = 2, respectively. 
We will briefly explain these problems here and refer the reader about vari- 
ational inequalities to |Glo84| and for two phase membrane to |Wei98] . 

• One phase obstacle problem (m = 1). Consider the following energy 
functional 

(2.2) minE{u)= j -\Vu\^ + fudx, 

over the convex set = {u £ H^{il) : u > 0,u = cp > on dO,}. The 
minimizer of ( 2.2 ) satisfies the following Euler-Lagrange equation 

{Au = fX{u>o} in ^, 
u = (l) on dn, 

u = |Vn|=0 inf]\{u>0}. 

• Two phase membrane problem (m = 2). 

Let fi : — 7- M, i = 1,2, be non-negative Lipschitz continuous 
functions, where is a bounded open subset of M" with smooth 
boundary. Let 

K = {ve W^'^{n) :v-ge Wo'^(f^)}, 
where g changes the sign on the boundary. Consider the functional 

(2.4) I{v) = j Q|Vw|2 + /imax(u,0) - /2min(u,0)^ dx, 

which is convex, weakly lower semi-continuous and hence attains its 
infimum at some point u £ K. In the functional (2.4) set 

gi=g^, 92 = 9', 



where = max(ibw,0). Then the functional I{v) in ( 2.4) can be 
rewritten as 

(2.5) /(ni, u,) = + ^ + + dx, 

where minimization is over the set 

S = {{ui, U2) G {H^{VL))'^ : Ul • U2 = 0, Uj > Ui = §{ on dQ, i = 1, 2}. 

The Euler-Lagrange equation corresponding to the minimizer u is 
given by ( |Wei98j ) . which is called two-phase membrane problem. 

/2 g^ f Alt = fiX{u>o} - f2X{u<0} in ^, 

\^ u = g on dQ, 

where T{u) = d{x E : u{x) > 0} U d{x G Q : u{x) < 0} n 17 is 
called the free boundary. 
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2.2. Numerical approximation of Problem (A) 

In this section we present our numerical schemes. Our numerical method 



uses the properties in Corollary 2.5, It means that if m(x) = 1, x S Br, then 
our scheme solves Aui = fi locally. For all x such that m{x) = 2, we should 
solve 

A{ui -Uj) = fiX{u,>Q} - fjX{uj>0}- 

To explain our method, first let m = 2 then we have 
(2.7) A(ni - U2) = /iX{«i>o} - f2X{u2>0}- 



Equation (2.7) shows that A{ui — U2) is bounded and consequently by clas- 
sical results for elliptic PDE we get that U1—U2 G C^'" for a < 1. Therefore, 
on the free boundary we have 

Vui = —VU2- 

Now for a given uniform mesh on C M^, let Uk{xi, yj) for k = 1,2, denote 
the average of Uk for all neighbor points of (xj, yj); 

Uk{xi,yj) = ^[uk{xi^i,yj) +Uk{xi+i,yj) + Uk{xi,yj-i) + Uk{xi,yjj^i)]. 

We use the standard finite difference discretization for equation ( |2.7| ), by 
letting Ax = Ay = h, to arrive at 

1 _ 1 _ 

^2 ^[4^xi(xi,yj) - 4ni(xi,yj)] - -^[4n2(xi, y^-) - 4u2(xi,yj)] 

= flX{ui{xi,yj)>0} - f2X{u2{xi,yj)>0}- 



We can obtain ni(xj, yj) and U2{xi,yj) from (2.8) and we impose the follow- 
ing conditions 

ui{xi,yj) ■ U2ixi,yj) = and ui{xi,yj) > 0, U2{xi,yj) > 0. 

Then the iterative method for ui and U2 will be as follows. 

• Initialization: 

,,(0). . _ r (xi,%) G n°, 

"1 [^^,VJ)-<^ Mx,,y,) (x„%)G5a 

,,(0). . _ f {xi,yj) £ n°, 

2 02(x„y,) (x„y,)G5a 

Here by Q° we mean the interior of the domain Q. 

• Step k + l, k>0: 

We iterate over all interior points by setting 

i[^'^'^\xi,yj) = max(^J^^^^^^^ +u[^\xi,yj) -u^^^\xi,yj), o] 



u''2^'^'^\xi,yj) = max(^^^%^^-F4''^(xj,yj)-nf^(xi,yj), o] 



Note that if m = 1 then the above method can be modified. The con- 
vergence of method in this case is given in [Glo84j . The algorithm for an 
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arbitrary m is as follows. Suppose there is a grid on the domain then our 
method can be formulated as follows: 

• Initialization: For l = l,--- ,m, set 

• Step k + 1, k > 0: For / = !,••• , m, we iterate for all interior points 

-fih' 



(2.9) 



max 



P^i 



Remark 1. Note that this iterative method is slow since information propa- 
gates from the boundary into the domain. One interesting question is, how 
can the idea of multi- grid method be applied? 



Lemma 2.6. The iterative method (2.9) satisfies 



,(k) 



r'{xi,y,).u'g'^\x,,yj) = 0, 
for all k £N and q, I £ {1, 2, . . . , m} , where q ^ I. 



Proof. Observe that from (2.9) it follows that 

uf\xi,yj) > 0, 



for all /c G N and / G {1, 2, . . . , m}. Assume u^^\xi, yj) > then by (2.9) we 
have 



ul {xi,yj) 



+ u\ ^'{xi,yj) -^u^p^ ^\xi,yj). 



This shows that 



u 



P^i 



Thus 



k 1 f h"^ 

uv: ^\xi,yj) < uf'^\xi,yj) < + ^uf~'^\xuyj), 



and after rearranging above inequalities we arrive at 



(2.10) 



+ 4^ ^\xi,y^)-^uf i)(x„yj)<0. 

V+1 



In light of (2.10) and (2.9) we derive 



u[''\xi,yj) = max 
Thus 



+ uf-i)(x„y,) -J2uit-'\x„y,), I = 0. 



Lf\xi,yj) ■ uf\xi,yj) = 0. 



9 



A. Arakelyan& F. Bozorgnia 



□ 



In order to see the consistency of the method to our problem (1.1), we 



wiU consider the finite difference scheme of our method (2.9). The scheme 
apparently will be the following discrete nonlinear system : 



(2.11) 



{Xi,yj) G 9n, 



where / G {1, 2, . . . , m}. 



Now we are going to show the consistency of the scheme (2.11 ) to the prop- 



erties that discussed in the Corollary |2.5[ First of all the disjoint property 



of the components follows directly from Lemma 2.6 Suppose ui{xi,yj) > 0, 
and together with this Up{xi, yj) = 0, for all p ^ I. This will imply that 



ui{xi,yj) 

and hence 
(2.12) 



-fih' 



+ ui{xi,yj) - '^Up{xi,yj) 



1 



{Aui{xi,yj) - Aui{xi,yj)) = fi{xi,yj). 



But equation (2.12) is just a discrete scheme of the Poisson equation 

Aui = fi. 

Hence if in the discrete sense u/(x, y) > 0, then we have Aui = fi- If we are 
locally on the free boundary of two components, say ui and Uq, then in the 



scheme (2.11) we have the following situation: 



ui{xi,yj) = Uq{xi,yj) = Up{xi,yj) = 0, where p / / and p q. 



According to the scheme (2.11) we have 



and 



= max ui{xi, yj) - Uq{xi, yj) 



= max (^Uq{xi,yj) - ui{xi,yj) 



4 



, 



Therefore 

4 _ 

-fq{xi,yj) < j^{ui{xi,yj) -Uq{xi,yj)) < fi{xi,yj), 
and by taking into account that ui{xi,yj) = Uq{xi,yj) = 0, we derive 

-fq < Ah{ui - Uq) < fi, 

at {xi,yj). Combining all results we will see the consistency with Corollary 
[231 
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Here we give a proof of the convergence of our method to the discretized 
problem, in the case m = 2 and fi = 0. We consider the following non-linear 
finite difference method 



(2.13) 



u 



k+l 



(xi, Vj) = max(tif - ii2, 0), 
U2^^{xi, yj) = max(ti2 — , 0). 



Note that (2.13) can be written as 
(2.14) 



h 

k+l. 



(xi, Vj) = max(uj - n|, 0) = ^ ( 



n'i 



U2 + 1^1 



k\ 



I Li2 \xi, yj) = max(u2 — u\, 0) = ^ — u\ + \u2 — u\\) . 
By subtracting the first equation from the second we will obtain 



(2.15) 



u 



Xi ) Vj ) 



Xi ) yj ) 



U2, 



which is a classical finite difference scheme for A{ui{x) — U2{x)) = 0. It 
is noteworthy that (2.15) follows from the last part in Corollary 2.5 This 
means that we have convergence of 

(2.16) u\{xi,yj) - U2{xi, yj) 



at every point {xi,yj), when k — oo. Recalling that by lemma (2.6) 
u\{xi,yj) ■ U2{xi, yj) = 0, for every > 0, we can write the following identity 
for all k. 



(2.17) 



{ui{xi,yj) -U2{xi,yj))'^ = {u'l{xi,yj) + U2{xi,yj))'' 



Therefore convergence of Ui{xi,yj) — U2{xi,yj) at every point {xi,yj) will 
imply the convergence of 

{ui{xi,yj) - U2{xi,yj))^, 



at every point as well. Hence, by (2.17) the sequence 

{ui{xi,yj) + U2ixi,yj))'^ 

converges at every point {xi,yj). Note that u\{xi^yj) and u\{xi,yj) are 
positive, this will imply the convergence of 

u\{xi,yj) + U2{xi,yj). 

Now convergence of u\{xi, y^)— u^Xj, yj) and Ui{xi, yj)+U2{xi,yj) will imply 
the convergence of Ui{xi,yj) and uKxj, ?/j) at every nodal point {xi,yj). This 
completes the proof. 



3. Theoretical results of Problem (B) 



In this section we present some result that have been proved for Problem 
(B). Here we mention some results for the case of two-species in dimension 
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two. Consider the following system. 



(3.1) 



ut — diAu = Au(l — u) 
vt — d2^v = Xv{l — v) 
u{x,y,t) = <t){x,y,t) 
v{x,y,t) = ^{x,y,t) 
uix,y,0) = uo{x,y) 
, v{x,y,0) = vo{x,y) 



1 2 
- ^UV"^ 

1 2 



in n X (0, oo), 
in X (0, oo), 
on dn X (0, oo), 
on dQ X (0, oo), 
in fi, 
in 17. 



This problems has been studied in !CDHn7b[ IEY941 ETM, |Squ08| , where 
the references of some physical background involving cubic coupling is given. 
This system, for steady boundary data admits a Lyapunov energy. Assume 
that the initial conditions no(x,y) and VQ{x,y) have disjoint supports and 

< uo{x,y),vo{x,y) < 1. 

We also assume that the boundary conditions are positive with disjoint 
support. The following Theorem has been proved in |Squ08] . 

Theorem 3.1. There exist two functions u{x,y),v{x,y) G H^{0.) n L^{Q) 
such that 

{ue^{tm),Ve^{tm)) {u,v) in ^{9) X ^{Vt) for any p > 2, 
as e — )■ and t — t- oo, where < u,v < I and u ■ v = in fl. Moreover, 



u 



-diAu < Xu{l 
and u[gn= </>, v[on= ^p. 

Now consider the following system 



-d2Av < Xv{l - v), 



(3.2) 



Ut — diAu = Xu{l — u) 
Vt — d2A.v = Xv{l — v) 
u{x,y,t) = (l){x,y,t) 
v{x,y,t) = ip{x,y,t) 
u{x,y,0) = uo{x,y) 
^ v{x,y,0) = vo{x,y) 



-uv 
1 

-uv 



in $7 X (0, oo), 
in $7 X (0, oo), 
on (90 X (0, oo), 
on (90 X (0, oo), 
in O, 
in O. 



It has been shown in j CDH+04j that for any T > as e tends to zero 
there exists a sequence of solutions [usiVe) to the system (3.2) converging 
in L^(0 X (0, T)) to a bounded segregated state (u, v), such that w = u — v 
solves the limiting free boundary problem (3.3), which shows the spatial 
segregation phenomena on finite time intervals. 

Theorem 3.2. iCDH+04j LetT > 0. Then there exists a sequence Em and 
u,v ^ L°° with 



{ue^.vej ^ {u,v) m L\n X (o,r)) xL\nx (o,r)), 
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as £ — )• 0, where < u,v < 1 and u - v = in Moreover, w = u — v is the 
unique weak solution to the following free boundary problem 

{wt — AD{w) = Xw{l — \w\) in Q X (0, oo), 

D{w{x, y,t)) = di(j){x, y, t) — d2i^{x, y, t) on dft x (0, oo), 
w{x,y,0) = uo{x,y) - vo{x,y) in Q, 

where 

(3.4) D{a) = 1 f'' ^ - °' 

^ ^ ^ ^ \ d2a a < 0. 

Also the cases of time-dependent boundary co nditions an d possibly dif- 
ferent diffusion coefficients has been discussed in |CDH+04| . In |CDH07a| 



in the case of equal diffusion coefficients di = d2 and stationary boundary 
conditions, Crooks, Dancer and Hilhorst studied the long-term segregation 
for large interactions. They reduced the system to a single parabolic equa- 
tion whose solution have e-independent uniform bounds. This system does 
not admit a natural Lyapunov functional and a direct analysis is therefore 
for long term behavior is not possible. 

3.1. Numerical approximation of Problem (B) 

We present a numerical scheme for elliptic system in Problem (B) as e — t- 
0. To explain the method, we assume that there exists two components. 



Theorem (3.1) shows 

d2l^v - diAu = An(l - n)x{„>o} - Aw(l - v)x{v>o}- 

This equation is solved numerically by employing second order, centered, 
finite differences on the given grid i.e, 

(3.5) - ^[4:u{xi,yj) -4:u{xi,yj)] + 'l^[4v{xi,yj) -4:v{xi,yj)] = 

\u{Xi,yj){l - u{xi,yj))x{u(x„y,)>0} - >^v{Xi,yj){l - v{Xi,yj))x{v{x„y,)>0}- 



It is easy to see that the equation (3.5) is a quadratic equation with respect 



to u{xi, yj) and v{xi, yj). By similar description of Section 2.2, if u{xi, yj) > 
then we set v{xi,yj) = and vice versa. Let 4a = Xh^ then from (3.5) we 
get the following iterative formulas 

.(^+^)(x„,,) = maxf 2(d,^(^)(..,y,)-d,iyW (..,,,)) 

di-a+ yj{di - a)'^ + Aa{diu^''\x^,y.j) - d2V^^\xi,yj)) 

2{d2V^^Hxi,yj) - diu''^\x^,yj)) 




d2-a + ^J{d2 -a)2 + Aa{d2v'-''^Xi,yj) - diu^^\xi,yj)) 



,0 



This approach can be extended for m components as well. The idea is 
just to take the difference between the i-th equation of the system and the 
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sum of all other equations. After that we use the same disjointness approach 
by setting Ui{xs,yr) > and Uj{xs,yr) = for all i ^ j, on the grid point 
{xs, Vr)- This will lead us to the quadratic equation w.r.t Ui{xs, y,.) as above. 
Thus according to the same arguments as above for m components we obtain 
the following iterative method: For alH = 1, . . . , m, 

(3.6) ur'\^.y,) = max ( ^^^''^^^^m) _ , , 

\di - a + y (d/ - a)2 + Aawi^''\xi,yj) 

where 

Wi~''\xi,yj) = diuf-'^\xi,yj) -^^dpU^^^\xi,yj). 



Again by using the same approach as in Lemma 2.6, one can prove the same 
result for this method as well. 



Lemma 3.3. If mind; > a, then the iterative method (3.6) satisfies 



uf\xi, yj) • uf\xi, yj) = 0, 
for all k and q, I G {1, 2, . . . , m}, where q ^ I. 

3.2. Parabolic case 



In the case that the coupling term is ™, Theorem (|2.6|) states that w = u — v 



solves the limiting free boundary problem in Theorem (3.3) which shows the 



spatial segregation phenomena on finite time intervals. To solve the problem 



(3.5) the second-order, implicit, Crank-Nicolson method is applied. 



w 



n+l 



{xi.yj) -w'^{xi,yj) 



-(ADu;|"+^ , + l^Dw\l. „ 

(^3 7) dt 2^ '(^"%) 

= ^['u;"+^(l - w'^+^) + 'u;"(l - u-")]. 

In this case we can obtain an iterative formula for i(;"^^(xj, yj) as a func- 
tion of w^{xi,yj), w^{xi,yj) and w^~^^{xi,yj). 

4. Numerical Examples 

In this section we present different examples for Problem (A) and Problem 
(B). We consider the following minimization problem 

(4.1) / = / y (hvuil^ + fiu^dx, 



i=i ^ 



over the set 5 = {(ui, . . . , Um) G {H^{Q))"^ '■ Ui > 0,Ui ■ uj = 0, Uj = (/){ on 
Examples 1, 2 and 3 show the numerical approximations of Problem (A) for 
different values of m and different O. 
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Example 4.1. The Figure |4] shows the solution of (4.1) for the case that 



n = l,m = 2. We choose /i = 2 + sin{x),f2 = 1 + ar. The equation for 
til — ii2 is as follows: 



(4.2) 



{ui - U2)" = (2 + sinx)x{ui>o} - (1 + x'^)X{u2>o}^ x e [-2, 2] 
ui{-2) = 1, U2{2) = 1. 




Figure 1. The plot of ui + U2- 



Example 4.2. Consider problem (A) with m = 3 and /i = /2 = 1, /s = 
The free boundary is shown in Figure [2] The boundary value g is given by 

{A-x^ -2 < X < +2 k y = -2, 
4-y2 -2<y<+2&x = -2, 
^ -2 < a; < +2 & y = -X. 




(a) The free boundaries (b) ui + U2 + U3 



Figure 2. The left picture shows the free boundaries of solutions. The 
right picture shows the surface of ui + U2 + U3. 
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(a) contour (b) mi + U2 + «3 + U4 



Figure 3. The left picture shows the contours of solutions and zero 
set. The right shows the surface of ui + U2 + + U4,. 



Example 4.3. Let Q = [-1,1] x [-1,1] and m = 4, /i = 8;/2 = 6; /s 
2; /4 = 1. Boundary values (pi, (i=l,2,3,4) are given as follows. 



l-3;2 X G [-1,1] & y = 1, 
elsewhere. 



l-y2 ye [-1,1] &x = l, 
elsewhere. 



4>3 



l-x"^ X G [-1,1] & y = -1, 
elsewhere. 



elsewhere. 



Example 4.4. Let be as in previous example and m = 4, /i = 8, /2 = 

— y^|,/3 = 0, /4 = |x + y|. The boundary condition cpi, (1=1,2,3,4) are 



the same as in example 4.3 The interfaces are shown in Figure 4. 





u 


1 




U2 






U4 












U3 





-1 -0.8 -D.6 -0.4 -0.2 0.2 0.4 0.6 0.8 



Figure 4. The picture shows the free boundaries of solutions of 
Ml, U2, W3 and W4. 
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Boundary values 




X 



Figure 5. Boundary values for u{x,y); left v{x,y); in the right. 

Now consider the following system of m differential equations for i = 
1, • • • , m, as e — )■ 0, 

{-diAui = Xui{l - Ui) - \ui{x) Y^-^i ™ ^> 

> in 
Uiix) = 4>i{x) on 9r2. 



Example 4.5. Let 17 = [0,1] x [0, 1], m = 2, A = 1, di = 1.5, (i2 = 1. The 
steady boundary values for u{x,y,t),v{x,y,t) are defined by 



0.5- 2. 5x 0<x<0.2, i+N_i0.5-|x < x < 0.2 



8- 



°' = I O' ■ 0.2 <x<l, 1' i 0.8 < X < 1, 

0(O,y,t) =0.5, 0(l,y,t) = O, 

and 

V^(x,0,tj-| _i^5^ 0.2<x<l, \ -2 + 2.5X 0.8<x<l, 

V'(0,2/,t) = 0, ^(l,y,i) = 0.5, 
Figure [5] shows boundary values of u and v. In the Figure [6] the contours 
plot of solutions u{x,y) and v{x,y) are shown. 
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Figure 6. Contours of u and v. 
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